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Abstract 

The nonequilibrium ionization process in 
hydrogen-helium mixture behind a strong shock wave 
is studied numerically using the detailed ionization rate 
model developed recently by Park which accounts for 
emission and absorption of Lyman lines. The study 
finds that, once the avalanche ionization is started, the 
Lyman line is self-absorbed. The intensity variation of 
the radiation at 5145 A found by Leibowitz in a shock 
tube experiment can be numerically reproduced by 
assuming that ionization behind the shock wave prior 
to the onset of avalanche ionization is 1.3%. Because 
1.3% initial ionization is highly unlikely, Leibowitz’s 
experimental data is deemed questionable. By varying 
the initial electron density value in the calculation, the 
calculated ionization equilibration time is shown to 
increase approximately as inverse square-root of the 
initial electron density value. The true ionization 
equilibration time is most likely much longer than the 
value found by Leibowitz. 

1. Introduction 

In the past, on-going, and future entry flight missions 
to the outer planets, heating rates to the vehicle’s 
surface tend to be large, and therefore the extent of 
ablation becomes also large. Therefore, accurate 
prediction of the heating rates becomes imperative. 
Computational Fluid Dynamics (CFD) can be a helpful 
tool in predicting the heat transfer rate if its reliability 
can be validated against a flight data. 

The atmospheres of outer planets consist of hydrogen 
and helium. An entry flight into the planet Jupiter has 
already been accomplished in the Galileo Probe 
mission. Surface recession data have successfully been 
obtained in that mission [1], This data is the only flight 
data available to validate the CFD methodology in 


designing the heatshield. 

The Galileo Probe data on surface recession has shown 
a surprisingly low surface recession in the stagnation 
region, and a surprisingly large recession in the 
downstream frustum region compared with the 
pre-flight predictions [2]. Very recently, Matsuyama et 
al [3] were able to explain the high surface recession in 
the downstream region by injection-induced turbulence. 
Flowever, the low surface recession in the stagnation 
region has not yet been explained. 

Park [4] has earlier speculated that the low heating in 
the stagnation region might be due to thermochemical 
nonequilibrium. The pre-flight predictions were made 
assuming thermochemical equilibrium. The 
nonequilibrium idea was first introduced by Howe [5], 
who derived his concept from the shock tube data of 
Leibowitz [6]. Leibowitz’s data showed that the region 
immediately behind the shock wave was not ionized 
and did not radiate. Howe derived an empirical 
expression for the thickness of the non-radiating 
nonequilibrium region, i.e. ionization equilibration 
time, in a hydrogen-helium mixture. 

In a companion paper to the present paper, Park [7] 
examined the nonequilibrium issue still further. He 
points out that the experiment by Leibowitz and 
subsequent experiment by his colleagues Livingston 
and Poon [8], which extended Leibowitz’s data to 
higher flight speeds, are mostly likely erroneous 
because of the interaction of radiation emanating from 
the arc-heated driver gas. Park [7] points out that the 
experiment by Stalker [9] in a shock tunnel produced 
an equilibration time 8 times longer than that by the 
Leibowitz-Livingston-Poon group. Park speculates that 
the true equilibration time may be even longer than 
that determined by Stalker. 

According to the experiment conducted by Bogdanoff 
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and Park [10], the freestream flow in front of the 
normal shock wave in a hydrogen-helium mixture is 
ionized by photo-ionization of H 2 . The temperature of 
the electrons so-produced is speculated to be 
determined indirectly by the temperature of the 
vibrational-rotational mode of H 2 through the 
electron-vibration-rotation coupling [7], The relaxation 
of vibrational-rotational mode of H 2 has been studied 
by Furudate et al [11], The study shows that the two 
modes are strongly coupled, and the temperature of the 
combined vib-rotational mode approaches the 
translational temperature relatively slowly, i.e., with a 
collision number of the order of hundred. The slowly 
rising electron temperature leads to slow initiation of 
the so-called avalanche ionization which sets off the 
rapid electron-impact ionization process that brings 
about the ionization equilibrium. The rate of avalanche 
ionization is dictated partly by the absorption of 
Lyman radiation [12]. Thus, the time for equilibration 
is dictated by the photo-ionization rate, vib-rotational 
relaxation rate, the coupling rate between electron 
temperature and vib-rotation temperature, and the rate 
of Lyman line absorption. 

It goes without saying that further experiments should 
be carried out to determine the true value of the time 
for ionization equilibrium in hydrogen. Along with 
such an experiment, a CFD modeling of the 
phenomenon would be required. The purpose of the 
present paper is to carry out such a CFD modeling, and 
investigate how the photo-ionization, vib-rotation 
relaxation, electron-vib-rotational coupling, and 
nonequilibrium ionization proceed in a hypersonic 
FL-Fle flow. Because of the enormous complexity and 
because most of the relevant parameters are unknown, 
a very first simple calculation is carried out in the 
present work. Flere, a new set of the ionization rate 
coefficients by Park [12] is employed. The ionization 
rate values used by Leibowitz in explaining his 
experimental data were arbitrary, and contained little 
physical ground. In comparison. Park’s new ionization 
rate values are firmly based on the state-of-the-art 
knowledge of such a process. The electron density and 
electron density were first chosen arbitrarily in order to 
numerically reproduce Leibowitz ’s shock tube data [6]. 
The parameters controlling the initial electron density 
and temperature were varied to show that the 
ionization equilibration times can be longer than that 
determined by Leibowitz. 


2. Methods of Calculation 


2.1 Governing equations 


The governing equations are the one 
dimensional Euler equations, 

dt dx V ’ 

In the equations, the global mass, the momentum, the 
total energy, the species mass, and the electron energy 
conservation equations are included. Hence, the 
conservative variables Q , the convective flux vector 
F , and the source term W are given respectively as 
follows; 
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where s stands for the chemical species, H 2 H , 
H * , He , He * , e . Electron-electronic energy E e , in 
the present study is defined by 

E el =1.5N e kT e +^E H (i)N H (i) ( 3 ) 

where N e is electron number density, k is 
Boltzmann constant, T e is electron temperature, E(i) 
is energy level of state i . Number density of atomic 
hydrogen in the state i , N„ ( i ) , is given by 


N H (i)= g-exp i-EJkT') 
expt-E, / kT e ) 


(4) 


In the present study, first three electronic states are 
considered. The electron-electronic energy source term, 
W el , can be written by 

ffi 3 

w e , = 2 N e V «,* — ——k(T — T e )—DW H 

tte m t 2 ( 5 ) 

+ X E H (0 N H (0 

where D is ionization energy and W H is chemical 
source term for atomic hydrogen. Elastic collision 
frequency, «*, is defined by °° ak = N k Q et -J8kT e / &n e . 
In the present study, the elastic collision cross section, 
Q et , are given by the same formula as in [6]. 


2.2 Chemical reaction rates coefficients 


The dissociation of molecular hydrogen, the ionization 
of atomic hydrogen, and the ionization of helium are 
considered. They are summarized in Table 1 . 
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2.2.1 Rates coefficients used by Leibowitz 

In [6], Leibowitz employs the idea of the two-step 
excitation-ionization process in his calculations; The 
atomic hydrogen is firstly excited by collisions 
(Reaction 3 and 5), and then ionizes rapidly [6]. The 
Electron-electronic energy Ed is defined by 
E d = 1.5 kN e T e in [6]. The reaction rate coefficients for 
these reactions are summarized in Table 2. The 
backward reaction rates are given by the fraction of the 
forward reaction rates and the equilibrium constants. 
The equilibrium constants are summarized in Table 3. 
Hereafter, this approach is called the Leibowitz ’s 
method. 


Table 1. Reaction 

Reaction 

1 

H + e °° H + +e + e 


2 

He+e 00 He + +e + e 


3 

H + e °° H‘+e 


4 

He+e °° He’+e 


5 

H + H “ H‘+H 


6 

H + He °° H‘+He 


7 

H 2 + He ~ H + H + He 


8 

H 2 + H 2 °° h + h + h 2 


9 

h 2 + h ~ h + H + H 


10 

H 2 + H* “ h + H + H + 


11 

H 2 +e o» H + H + e 


Table 2. Reaction rate coefficients [6], 

Reaction 

Reaction Rates [m 3 /mole-s] 

Ref. 

1 

6.09 =1 0" 3 y]8kT e / exp(- 15782/7;) 

13 

2 

3 . 56 ~i o- 23 ^/8 kT e /*°c exp(-2 85248 / T e ) 

14 

3 

7.5 ~1 0" 2 j8kT e / exp(- 1 1 605 / T e ) 

15, 16 

4 

6.00=1 0" 3 V8 kT„ T#> 7 exp(-23210 IT e ) 

17 

5 

4.0=1 0 -23 78AT / exp(- 1 1605/T) 


6 

4.0=10" 3 V8AT/^=exp(-11605/r) 


7 

6.93=10 _12 /T exp(-52340/T) 

18 

8 

2.5k, 

18 

9 

20.0 k, 

18 

10 

20.0k, 


11 

20.0k, 


Table 3 Equilibrium constants [6]. 

i Equilibrium constants, K, [mole/nr 1 ] 

H + 

4.05=10- 3 r 3 ' 2 exp(-15872/T) 


He * 

1.62~10- 2 r 3 ' 2 exp( -285248 /T) 


H 

3.7~10 <i [1.0-exp(1.50=10 8 r _2 )]exp(-52340/r) 


2.2.3 Park’s ionization rate coefficients of atomic 
hydrogen accounted for radiation absorption 

In order to account for the effect of the radiation 
absorption, the rate coefficient for the electron-impact 
ionization (Reaction 1) is taken from [12], and 
Reaction 3 is discarded. The change of electron 
number density due to the electron impact ionization of 
atomic hydrogen is expressed by 
dN 

— L = k f N H N e — Nl°° x — /V/°?“fl) , (6) 

dt 

<*>! = =— “(1) = N e k r + V =(;) , (4) 

i=2 

where k f is the collisional ionization rate coefficient, 
k r is the collisional three-body recombination rate 
coefficient, N H is the atomic hydrogen number 
density, = is the collisional-radiative recombination 
rate coefficient, «(/) is the rate coefficient for 
radiative recombination into state i , <* c is Lyman 
continuum radiation fraction, and m is the highest 
bound state quantum number. Using the source code 
provided in [12], the parameters, k , , k r , = and 
oc(i) can be obtained from the number density N H 
and N e , the electronic temperature T e , and the Lyman 
line absorption factor , The Lyman °° line 
absorption factor is defined by 

_ 5(1,2) _ P/1.634~10- 11 iV / ,(l) 

~ 5(1, 2) e ~ 2 2 exp[- E h (2) / kT c ]A(2, 1) ' ^ 

where 5(1, j) is rate coefficient for radiative transition 
from lower state i to upper state j and subscript E 
stands for a equilibrium state. Power absorbed by the 
Lyman <= line, P , can be written by 

P= 1.634=10-“ 4(2, l)iV(2) - °°q md , (6) 

where °°q rad is divergence of radiative heat flux 
determined by solving the radiative heat transfer 
equation. 

2.2.4. Initial electron density and temperature 

Initial electron density and temperature at the onset of 
avalanche ionization were varied arbitrarily because 
they are totally unknown at this time. The relative 
initial electron concentration (N e /N H2 ) 0 was varied 
between 5=10^ to 0.05 , and the initial electron 
temperature was varied from 300 to 6000 K. The true 
initial electron density value is most likely lower than 
the lowest value considered here. However, lower 
initial electron density values led to numerical 
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instability. This problem needs to be solved in the 
future. 

2.3 Numerical method 

The convective fluxes are given by AUSM-DV upwind 
scheme. The dependent variables are interpolated by 
the MUSCL approach An explicit Runge-Kutta method 
is used for time integration. 

3. Test condition 

A shock tube problem is numerically solved for the 
experimental condition of Leibowitz. 3 The test gas is a 
mixture of 20.8% H 2 -79.2% He. The driver gas is He. 
The measured shock velocity was 15.5 km. Using this 
value of shock velocity, it was impossible to 
numerically reproduce the measured electron density: 
the measured electron density values were higher. The 
measured electron density values were numerically 
reproduced assuming the shock velocity to be 16.6 
km/s. The reason for such a high effective shock speed 
is believed to be the radiative heating of the driven gas 
by the radiation emitted by the driver gas [10]. The 
static pressure at the test section is 133 Pa. The initial 
translational and the electron temperature in the test 
section are set to 300K. The pressure ratio between the 
test section and the driver section is determined from 
the Rankine-Hugoniot relation. The computational 
space is 2.0 m in length and divided by 1000 cells with 
a constant space. The first 120 cells are considered to 
be the driver section. 

4. Results and discussions 
4.1. Leibowitz data 

In order to verify the present numerical code, the 
reproduction of the calculated results by Leibowitz is 
attempted using the same reaction rates as in 
Leibowitz ’s [6], The initial electron density is assumed 
to be 0%: avalanche ionization is assumed to be 
initiated by electrons produced mostly by Reaction 5, 
as in Leibowitz’s work. In Fig. 1, the calculated 
temperature profiles behind the shock wave are 
compared with the Leibowitz’s numerical result. In the 
present calculation, the position of shock wave is 
defined by the position where translational temperature 
begins to increase. Differences in the temperature 
variation between these two calculations are initiated 
immediately behind the shock wave. They may be 


caused by differences in the numerical methods. Note 
that, in Leibowitz’s calculation, the equations of fluid 
mechanics are solved in a shock fixed coordinate 
system, and convections of fluids are ignored. A good 
agreement in the temperature is obtained in the 
downstream where the transnational and the electron 
temperature are in equilibrium. Fig. 2 shows the 
corresponding number density profiles for atomic 
hydrogen and electrons. Deviation in the atomic 
hydrogen number density is also believed to come 
from the differences in the numerical methods. 

In order to qualitatively examine the effect of radiation 
absorption on ionization of atomic hydrogen, 
calculations with Park’s ionization rate at three fixed 
values of the Lyman line absorption factor », are 
implemented. The values of °v = 0, 1 and 100 are 
chosen, representing an optically thin case, a 



Fig. 1 Temperatures profile obtained in the present 
calculation. 



Fig. 2 Number density profile obtained in the present 
calculation. 
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Fig. 3 Effect of Lyman radiation on electron number density 
profde. 
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Fig. 5 Radiative power absorbed by the Lyman = line and 
Lyman line absorption factor behind shock wave. 




Fig. 4 Effect of Lyman radiation on temperatures profde. 


Laboratory time, [°<s] 



Fig. 6 Calculated radiative intensity and the experimental 
data by Leibowitz (==5145 A). 


self-absorbed case, and a strongly irradiated case, 
respectively. For all the cases here, the Lyman 
continuum radiation fraction «. is set to 1. Fig. 3 
shows the number density profiles obtained in the 
present calculations. The electron number density stays 
low in the first few microseconds behind shock wave, 
where the collisional ionization (Reaction 5) is 
dominant. Then, an avalanche ionization zone appears 
where the electron-impact ionization is dominant. In 
the cases for higher , incubation time to initiate the 
avalanche ionization is shorter, and the rate of electron 
production during the avalanche ionization is faster. 
The electron number density at the quasi-steady state 
(QSS), which appears after the avalanche ionization 
ends, is also affected by . In the strongly irradiated 
case of °y = 100, the electron number density at the 


QSS is 30 % larger than the case of <* L = 0. The 
temperature profiles are also drastically changed 
according to the value of <* L , as shown in Fig 4. The 
high degree of electron impact ionization with the large 
value of <* L lowers the electron temperature in the 
ionization zone. Accordingly, equilibrations between 
the transnational and the electron temperature proceed 
faster. 

The variation of the value of « c gives a miner effect 
on the flowfield. Although the figure is not shown here, 
the electron number density profiles with «. = 0.5 are 
almost identical to the results with <*. = 1 when the 
value of <* L is fixed. 

A calculation using local values of <* L , which is 
determined by solving the heat transfer equation, is 
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performed next. The heat transfer equation is solved 
using a multi-band model. Absorption coefficients for 
the Lyman °o line are evaluated at 23 wavelength points 
in the range from 1200 to 1230 A. The calculated 
radiation power of absorption by the Lyman °o line and 
the Lyman line absorption factor <* L arc shown in Fig. 
5. The radiation power of absorption varies from 0 in 
the non-ionization zone behind the shock wave to 
about 2000 MW/m 3 in the QSS zone in the down 
stream. The Lyman °° line absorption factor ^ keeps 
constant at the value nearly equal to 1 after slight 
increase immediately behind the shock wave; the 
emission from Lyman °o line is totally absorbed. 
Therefore, the electron number density and the 
temperatures profiles are almost identical to those in 
the case of = I in Fig 3 and 4. The calculated 
radiative intensity profile for the wavelength of 5145 A 
is compared with the experimental data in Fig. 6. The 
same approach as in [3] is employed to calculate the 
intensity. The obtained intensity profile agrees well 
with the experimental data measured by Leibowitz. 

4.2. Effect of initial electron density and 
temperature 

In this subsection, effects of initial electron number 
density and electron temperature are examined by a 
parametric study. For the calculations in this study, 
Reaction 5 and 6 in Table 1 is discarded, and the initial 
number density at the point of initiation of avalanche 
ionization is given at a certain value. From the 
observation in the previous subsection, ^ = 1 and 
= 1 are employed for all the cases in this 
subsection. 

First, the ratio of the initial electron number density to 
the initial molecular hydrogen is varied, while the 
initial electron temperature is remained to be 300K. As 
shown in Fig. 7, the intensity profile for 5145 A is 
changed according to variation of (N e /N H2 )„ . 
Leibowitz ’s experimental data can be reproduced with 
( N e /N H1 ) 0 = 0.013. The initiation and the degree of 
avalanche ionization of atomic hydrogen are greatly 
affected by the initial electron number density, as 
shown in the electron number density profile in Fig. 8. 
Larger non-ionization region is obtained with lower 
initial electron number density. 

Next, the initial electron temperature, To , is varied 
from 300 to 6000 K, while the initial (N e / <V„ 2 ) 0 is 


Laboratory time, [<*s] 



Fig. 7 Effect of initial electron number density on intensity 
profile. 



Fig. 8 Effect of initial electron number density on electron 
number density profile. (To = 300 K). 



Fig. 9 Effect of initial electron temperature on electron number 
density profile ( (N c / N hi)u =0.013). 
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fixed at the value of 0.013. The initial electron 
temperature at the point of initiation of avalanche 
ionization is at presently totally unknown. The electron 
number density profiles obtained for the various value 
of Teo are shown in Fig. 9. For 2000 <T e o< 5000K, 
times to initiate the avalanche ionization are almost 
identical, although higher temperature gives more 
production of electrons. The time to reach QSS with 
T e o= 2000 is about two times large than that of 300 K. 
For the T e o higher than 2000K, the time to QSS 
becomes shorter as the temperature increases. It is 
interesting to see that the effect of high initial electron 
temperature is not monotonic: the longest delay of 
ionization occurs at electron temperature of 2000 K. 

Fig. 10 summarizes the effect of the initial electron 
number density on the e-folding ionization equilibrium 
time. For all the initial electron temperature cases, the 
equilibration time increases according to the 
(N e / iV,, 2 )„ decrease. When the initial electron 
number density is lower, the effect of initial electron 
temperature is miner. The equilibration time for 
(N e / N H1 ) 0 = 0.0005 is about 5 times longer than that 
for {N e I N H i)»— 0.013 at 7%=300K. For To =3 00 K, 
the ionization equilibration time varies approximately 
as inverse square-root of the initial electron density. 

5. Discussion 

In this study, only the avalanche ionization process is 
modeled. We find first that, once avalanche ionization 
process is started, the absorption of Lyman lines is of 
little consequence, because it is nearly totally 
self-absorbed. However, we do not yet know how 
Lyman radiation will affect the region prior to the 
avalanche ionization. 

Another important finding here is that the time to reach 
ionization equilibrium is affected strongly by the initial 
electron density and temperature at the point of 
avalanche ionization. In the range of the initial electron 
density values considered in the present work, a factor 
of 5 variation is seen. As mentioned in 2.2.4, the 
lowest value of initial electron density was set in the 
present work by the numerical instability. The true 
initial electron density is most likely lower than the 
value considered in the present work. In the region 
ahead of avalanche ionization, photo-ionization 
process is the only mechanism to produce the initial 
electrons. Because electron temperature will be low in 



Fig. 10 Effect of initial electron number density and initial 
electron temperature on the ionization equilibrium time. 


the region, the free electrons produced by 
photo-ionization will be rapidly recombining. The 
present work finds that Leibo witz’s shock tube data 
can be numerically reproduced with an initial 
(N e / N H2 ) 0 value of 0.013, i.e. 1.3% ionization. This 
value of 1.3% is much higher than the value found 
experimentally by Bogdanoff and Park [10], and 
should be considered impossible. The only explanation 
is the photo-ionization by the strong radiation 
transmitted from the arc-heated driver gas [10], 

According to Fig. 10, equilibration time increases 
roughly as an inverse square-root of the initial electron 
density, as mentioned. If the true value of (N e / N H2 ) 0 
is 10~ 6 , the equilibration time will be two orders of 
magnitude longer than Leibowitz’s value. The 
equilibration time value of Stalker [9], which is 8 times 
larger than the value by Leibowitz, is well within the 
explainable range. This means that, at this time, we do 
not know how to predict the thickness of the 
nonequilibrium, un-ionized, non-radiating region 
behind the shock wave in outer planet entry flights. 
This points to the need to know more about the 
processes occurring in the region prior to the initiation 
of avalanche ionization. 

6. Conclusions 

Using the detailed ionization rate model developed 
recently by Park, it is found that the absorption of 
Lyman alpha has a great influence on the initiation and 
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the rate of avalanche ionization behind a shock wave. 
For a strong absorption, avalanche ionization initiates 
earlier and faster. Once the avalanche ionization is 
started, the Lyman line absorption behind shock wave 
is almost constant at the value of 1, which represents 
the emission from Lyman line is self-absorbed by the 
atomic hydrogen gas. The calculated intensity profile 
for 5145 A agrees well with the experimental data by 
Leibowitz, when the starting degree of ionization is 
assumed to be 1.3%, which is unrealistically high. The 
ionization equilibration time is roughly inversely 
proportional to the square-root of the initial electron 
density, which are presently unknown. The true 
ionization equilibration time could be up to two orders 
of magnitude longer than that determined by 
Leibowitz. 
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